Simulating nanoscale dielectric response 
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We introduce a constrained energy functional to describe dielectric response. We demonstrate that 
the local functional is a generalization of the long ranged Marcus energy. Our re-formulation is used 
to implement a cluster Monte Carlo algorithm for the simulation of dielectric media. The algorithm 
avoids solving the Poisson equation and remains efficient in the presence of spatial heterogeneity, 
nonlinearity and scale dependent dielectric properties. 
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The response of dielectric materials leads to important 
modifications of the bare electrostatic energy of charges 
P|. Even on microscopic scales many important quali- 
tative insights result from the application of the macro- 
scopic laws of electrostatics. For example, a large part of 
the solvation free energy of an ion in water is understood 
in terms of the Born self energy of an ion, q'^ /Snae where e 
is the macroscopic dielectric constant of the solvent and a 
an atomic scale; similarly many essential features of ion 
channels can be explained by a continuum description 
of high dielectric channels through low dielectric mem- 
branes [113. 

On the nm-scale the dielectric response is in general 
non-local Dr = / er,r'Er' (fir' , where D is the electric 
displacement and E the electric field. In Fourier space 
the field energy of a homogeneous fluid is given by 

where e{q) is now a scale dependent dielectric constant 
with D(g) = e{q)'E{q). If we take water as an example 
we learn that starting from e ^ 80 at q = 0, the dielectric 
constant of water drops to e = 20 at a wavelength which 
is comparable to the Bjerrum length of a monovalent ion 
(7 A); below a wavelength of ~5A e diverges and be- 
comes negative. Moreover, many situations of practical 
interest require the inclusion of non-linear effects such as 
dielectric saturation or surface ordering and alignment 

Recently, we introduced an algorithm which permits 
the calculation of electrostatic interactions in heteroge- 
neous JocaJ dielectric media without solving the Poisson 
equation 0. The idea is to generate long ranged elec- 
trostatic interactions dynamically via local interactions 
between charges, the medium and the electric field. On 
first sight, it seems straightforward to replace the field en- 
ergy U = /Dr/2er d^r used in 0] by eq. iQ}. However, 
this approach is fundamentally flawed: the method would 
allow one to treat unphysical systems with < e(g) < 1 



while becoming unstable for real materials in which 
the dielectric constant becomes negative. 

In this letter we present a generalization of the Marcus 
energy for polarizable media Q written in terms of the 
truCjindependent thermodynamic variables in the prob- 
lem |j,|9|: the electric polarization of the medium, P, and 
the displacement, D. We show that our formalism is ca- 
pable of reproducing the full range of the linear dielectric 
response seen in nature including negative dielectric con- 
stants. 0, . In addition, we demonstrate that the same 
techniques can be used to study non-linear or polar di- 
electrics at essentially identical computational cost. The 
key feature is the use of a cluster algorithm for the 
equilibration of the field degrees of freedom. Similarly 
to related problems for classical ^3 ^^'^ quantum spins 
|l2l | the numerical effort per autocorrelation time can be 
reduced from 0{L^) sweeps with z > to 0{L'^) sweeps 
The prefactor in this scaling depends on the dielectric 
properties. 

We now introduce the energy functional. There are 
two contributions to the energy of a dielectric medium. 
Firstly, the energy density = (D-P)^/2 of the elec- 
tric field; secondly, an elastic polarization energy, G(P), 
due to short ranged interactions between molecules. We 
start by expressing G as a general quadratic function of P 
with a short ranged kernel Kj-y which we describe more 
fully below. In a heterogeneous system K varies from 
place to place in the sample; it contains all the material 
properties. Thus, 

We use units with eQ — 1 and periodic boundary con- 
ditions. Furthermore, D is constrained by Gauss' law, 
div D — p = 0. In the following we demonstrate the 
equivalence of our formalism to the standard theory of 
dielectric media [ll[ll[l3. 

We first work at zero temperature; this will allow us to 
calculate the relationship between the dielectric constant 
e{q) defined in eq. and the kernel K{q) in eq. We 
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minimize eq. ^ subject to the constraint of Gauss' law 
with the help of a Lagrange multiplier (p. We consider 
the stationary points of the functional 



A = U- 



j d^r (j){dWD-p) 



(3) 



We will pass rather freely between the full integral form, 
eq. (O and an operator notation in which all components 
of a field are grouped in a single vector and is a matrix. 
The variational equations are: 



SF 
5T> 



div D - p = 
P D + = 
D - P + grad(j!) = 



(4) 



We indeed find from the variational equations that if we 
define E = D — P then the relation between the polariza- 
tion and the electric field is P = xE, where x = is 
the susceptibility to the electric field, E = —gia,d(j>. We 
now solve eq. (@J for P in terms of D and substitute in 
eq. 121. We find eq. (QJ with 



eiq) = l + K~\q) 



(5) 



We have reproduced all the standard relations between E, 
P, (j) and D of electrostatics, as well as the energy eq. (Q) 
which should be the minimum of the functional eq. 1^)1. 
Similarly, we can show that our Ansatz is equivalent to 
the Marcus polarization functional 0, 0, Q|: we use 
eq. to eliminate the constrained field D and express 
eq. (0 in terms of the polarization P and the bare electric 
field Eq = — grad</)o generated by the free charges in 
vacuum: V^oin 



-P- 



r — r' 



,3^ .3„. divP.divP,, _ / , p ^3, 

Pri^r.r'Pr' d^T d^v' + I ^ d^T 



(6) 



It is instructive to rewrite Eq. © in the operator nota- 
tion: 



Up = -P(r + if)P-Eo-p 



E2 



(7) 



T denotes \q){q\ —^^^j^^f-^ + 5^(r) in real space. The 
field P interacts with itself via the long ranged dipole- 
dipole potential. 

We now consider the stability of the functional eq. Q 
for g ^ 0, in order to check that the energy is a true min- 
imum, not just a stationary point. The longitudinal con- 
straint of eq. (jSJ freezes D; only fluctuations of P are free. 
The coefHcient of P^ in eq. is V,(P) = F{l + K)P/2. 
In order for fluctuations in P to be bounded (so that 
the ground state is stable) we require that eigenvalues of 



(1 + K) are positive. Now express Vg(P) in terms of the 
field E and eliminate K for e using eq. lO): 



K,(E) 



> 



(8) 



Stability implies that e{q){e{q) — 1) > 0, so that e > 1 or 
e < 0. The expression eq. ® is indeed that given in~j3 for 
the effective potential of the electric field. Stable modes 
with e{q) < (such as those seen in water) correspond 
to —1 < K{q) < 0. We conclude that our constrained 
energy eq. 101 is capable of reproducing the full range of 
dielectric response seen in nature and leads to the cor- 
rect band of forbidden response functions < e(g) < 1 
where the system becomes thermodynamically unstable: 
in electrostatics there can be no equivalent of diamag- 
netic response, which might seem plausible from the con- 
ventional energy eq l(T]l. 

We now generalize to finite temperatures and study 
the partition function 



Z = 



J dP dT> e-^^ W (5(div D(r) - p(r)) (9) 



We shall integrate over either P or D to find the finite 
temperature generalizations of eq. and eq. 0: The 
P degrees of freedom are Gaussian, when we integrate 
over them we find the constrained partition function 



dDe-'^^S^'^''-J|5(divD-p) (10) 



where we have dropped unimportant numerical prefac- 
tors. This constrained integral over D was studied in 
detail in |3| • When dielectric properties are local, so that 
K{y,y') — k{y)I6{y — r'), eq. 1)10(1 describes particles in- 
teracting through an electrostatic potential which is a so- 
lution to the Poisson equation, div (egrad(/)) — — p, with 
e(r) = 1 + 1/K(r). In addition it gives the Keesom po- 
tential between fluctuating, classical dipoles. 

Rather than integrating over P we integrate over D in 
eq. 0: the (5-function constraint is imposed using the 
identity 2'k5{x) = J e"^^ d0. The integral over D is then 
Gaussian, as is that subsequently performed over (j). We 
find Z = j' dPe~^^f . Up is the energy of eq. 0. This is 
our principal formal result. It shows that with the energy 
functional eq. Q integrating over the constrained field D 
generates results which are equivalent to using the long 
ranged Marcus functional eq. 0. 

The treatment of the stability of the field P requires 
generalization at finite temperature: We must distinguish 
between the longitudinal and transverse parts of the op- 
erator K when q^Q: K^f. From eq. (|2Jl we calculate the 
fluctuations of the polarization field. Longitudinal fluc- 
tuations give f3S{q) = P{5P ■ SP)i{q) = + Ki{q)) 
whereas for the two transverse modes (3{SP ■ (5P)t(q) = 
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2/Kt{c\), so that the eigenvalues of Kt must also be pos- 
itive. Combining these two expressions we find the fluc- 
tuations in P for small q: 



P{5V ■ &V) 



1 



(2e+l)(e-l) 



{l + Ki) Kt 
= 3(e-l) 



K{0) 



(11) 



where we have used the fact that K becomes isotropic 
at small q. These expressions are familiar from the Kirk- 
wood theory of dielectrics li'l and will be used to extract 
dielectric constants from our simulations. Using these 
expressions is always numerically challenging since they 
link the dielectric properties to fluctuations in the polar- 
ization field. Accurate results require simulations which 
are many times longer than the equilibration time. 




FIG. 1: Dielectric constant deduced from polarization fluc- 
tuations for the dielectric eq. 11311 . L = 64. Data plotted 
as a function of |q| with |qp — 2^^{1 — cosgi). {"i^ni) are 
-I- (5,-2.3), o (2.5,-1.5), (1, -1.0), A (0.5,-0.9). Four days 
simulation per curve on a Pentium-4 processor. 

By choosing appropriate kernels K we can introduce 
semi-microscopic models of dielectric media with a great 
variety of dielectric properties. A systematic approach is 
to use a Landau-Ginsbur g ex pansion of the free energy 
for the polarization field The lowest order terms 

(corresponding to a linear theory G — VKV /2) are given 
by 

G=\j {kP^ + Kt(curlP)2 + KKdivP)2} Sv (12) 

Eigenvalues of Ki are (k -f Kiq'^), transverse eigenvalues 
are [n + Ktq^). The dispersion laws spht at 0{q^). Stabil- 
ity implies that k, and Kt must all be positive. From 



eq. (O the dielectric constant is e(g) = {1 + 1 / {k + niq^)} . 
It falls off monotonically with wavevector. 

Until now we have worked with linear media, which 
as we have shown analytically reproduce the standard 
theory of dielectric media. However our approach is not 
limited to linear models. Consider soft Langevin dipoles 
where the length constraint is imposed by an energy func- 
tion G = J 7(P^ — Pq)^ <fr. Dipoles with long ranged 
interactions described by the operator T in eq. ^ have 
long been used 0| to describe polar solvents. In contrast 
to linear models, they exhibit a saturation of the polar- 
ization degrees of freedom when interacting with strongly 
charged objects. To generate, in addition, scale depen- 
dent dielectric effects we include derivative terms in the 
free energy: 



G = 



K;(divP)^ a(graddivP)^ 



(13) 

We are particularly interested in the case k; < in order 
to produce systems with e{q) < 0, together with a > 0, 
necessary for stability at large wavevectors. At interfaces 
other contributions to G such as / grad^ • P d'^r select a 
preferred direction for the polarization and can be used to 
study ordering by surfaces. We leave such considerations, 
however, to future work. In our first simulations we used 
a = -0.4k,, 7.5/7. 

We discretize the fields on a simple cubic lattice of side 
L so that E and P are associated with the 3L^ links. The 
constrained field D is sampled by a cluster (worm) algo- 
rithm [let Il2| invented to study quantum spin models. 
We sample P with the Metropolis algorithm. The worm 
modifies a large cluster of 0{L^) D variables while con- 
serving divD. We define a sweep as Monte Carlo 
tries for P and one worm for D. We measure PS{q) and 
determine e(q) = 1/(1 — l3S{q)). The results are plot- 
ted for several sets of parameters in Figure ^ Passage of 
PS{q) through 1 gives rise to poles in the dielectric prop- 
erties, leading to some of the major qualitative features 
known in water: Firstly, a long wavelength dielectric con- 
stant satisfying e(0) > 1, secondly, a band of wavevectors 
with negative dielectric constant, thirdly, convergence of 
e to 1 for large q. 

A central point of this letter is the demonstration of 
the computational efficiency of our approach. In the fol- 
lowing we compare a linear, non-local model eq. H12|l . 
a heterogeneous, linear, local model, and a non-linear 
model consisting of soft Langevin dipoles eq. (fT!^ We de- 
termine equilibration times with a blocking method |l9l | : 
Starting from N — 2" recordings, x{t), one calculates an 
estimate of the mean and standard error for blocks of 
data of length 6 = 2"* with < to < rt. We studied the 
dynamics of variable x{t) = P^(i, q — 0), used to measure 
the (7 = dielectric constant through eq. Hll|l . Block- 
ing analysis leads to a monotonically increasing estimate 
of the standard error, in (x), <j{b). When b the block 
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size has reached the autocorrelation time of the simula- 
tion the standard error converges to cr^ = 2Tint{5x'^) / N , 
where Tint is the integrated autocorrelation time. We cal- 
culated the blocking curves for various values of system 
size and model parameters, Figure [3 Systems of dif- 
ferent size L have blocking curves which superpose with 
no rescaling of the data; the relaxation time (in sweeps) 
of the algorithm is independent of the system size. The 
cluster algorithm is indeed characterized by a dynamic 
exponent z — Q and not hindered by the heterogeneous 
material properties. When we modify the dielectric prop- 
erties we need to rescale both horizontal and vertical axes 
to superpose data. We find that the scaling variables are 
Na^{b)/xo{Sx'^) as a function of b/xo where xo is the 
q = susceptibility; for the heterogeneous system we 
used the appropriate xo for the large e region. Thus the 
simulation is slowed by a factor xo when dielectric prop- 
erties change. For Xo = 1 Tint 2 sweeps when simulat- 
ing with a simple quadratic energy for the polarization 
fluctuations. When using the soft Langevin dipole the 
simulation is approximately six times slower. 

The appearance of a long time scale in high dielectric 
media can be understood as being due to the large ratio of 
the amplitude of transverse to longitudinal fluctuations 
2(1 -I- K~^) = 2e, which splits the experimental longi- 
tudinal and transverse relaxation time scales |2flj |. This 
same splitting limits the efficiency of our algorithm. 0(e) 
sweeps are needed to fully equilibrate the system and to 
generate the Keesom - van der Waals interactions. Note, 
however, that the autocorrelation time of the longitudi- 
nal modes which are important for interactions between 
charges remains 0(1) sweeps even when the transverse 
and q — modes are slow. 




FIG. 2: Scaled statistical error as a function of scaled blocking 
factor. Curves for L = 10, 13, 18 and simulations of A'' = 2^^ 
sweeps. Linear: eq. 11211 . k = 1.0, 0.1, 0.03, 0.015, Kt ~ 0, 
Hi = k/S:. Non-linear (blue): eq. 11311 . 7 =2, 1, 0.4, a = k; = 
scaled with xo- Heterogenous linear with L = 18 where 
K = 0.25 for < z < 1/2 and k = 0.02 for L/2 < z <L scaled 
with Xo = 50. 

We conclude that our approach allows the investigation 



of electrostatic interactions in heterogeneous, non-linear 
and non-local dielectric media. The small computational 
costs for treating more complicated (elastic) polarization 
energies are due to our local formulation of the problem. 
Long ranged interactions are not calculated explicitely, 
but generated through a constrained energy functional 
for the electric field. This is in marked contrast to the 
conventional global solution of the Poisson equation. The 
equivalent minimization of the functional eq. lO can be 
written as P = (T + K)~^'En). Since updating the in- 
verted operator (T + K) in computer simulations is very 
time consuming |l3], most implicit solvent schemes rely 
on approximate solutions of the macroscopic laws of elec- 
trostatics in heterogenous dielectric media ji^. Experi- 
ence [23, li^l shows that the inherent advantages of the 
formulation in terms of local fields are not lost in off- 
lattice Molecular Dynamics implementations. We there- 
fore believe that the present work opens the way to the 
development of more realistic implicit solvents models for 
(bio)molecular simulations. 
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